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Accomplishments: 

A.  Development  and  Application  of  Liouville-von  Neumann  Dynamics  Method  for 
Efficient  Direct  Dynamics  Calculations 

We  developed  a  novel  first  principles  Molecular  Dynamics  scheme,  called  LvNMD,  based 
on  the  Liouville-von  Neumann  equation  for  density  matrix  propagation  and  the  Magnus 
expansion  of  the  time  evolution  operator.  The  scheme  combines  formally  accurate  quantum 
propagation  of  electrons  represented  via  density  matrix  P(t): 


^  =[F(t),P(t )] 
dt 

and  a  classical  propagation  of  R  nuclear  degrees  of  freedom: 
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where  F  =  h  +  G(P)  is  the  Fock  operator,  M  is  mass  of  nuclei.  E  is  the  potential  energy  for 
nuclear  motion  obtained  from  electronic  structure  as  a  sum  1 -electron,  2-electron  and  nuclear 
terms: 


E  =  Tr[hP ]  +  ±7V[G(P)P]  +  VNN  (3) 

The  solution  to  Eq.  1  is  obtained  via  Baker-Campbell-Hausdorf  (BCH)  expansion. 

2  3 

P(t)  =  P(0)-it[F,Pm-t-[F,[F,Pm]  +  it-[F,[F,[F,Pm]]  +  -..  (4) 

Some  advantages  of  LvNMD  are:  (a)  algorithm  is  free  of 
constraint  and  fictitious  parameters  (different  from  ADMP 
(Atom-centered  Density  Matrix  Propagation)  and  Car- 
Parrinello  methods),  (b)  it  avoids  diagonalization  of  the  Fock 
operator,  and  (c)  it  can  be  used  in  the  case  of  fractional 
occupation.  The  last  point  is  very  important  and  it  allows  to 
effectively  incorporate  an  ensemble  of  many  states  of  nearly 
same  energy  (as  found  in  transition  metal  clusters)  via 
electronic  temperature  (Fermi-Dirac  weighted  density  matrix). 

The  algorithm  has  been  found  to  be  very  stable  and  have  a 
very  good  conservation  of  the  total  energy  (see  Figure  1  for 
comparison  of  energy  conservation  simulation  of  our  LvNMD 
with  BOMD  (Born-Oppenheimer  MD)  and  ADMP  (Atom- 
centered  Density  Matrix  Propagation)  methods  for  a  complex 
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Figure  1 .  Conservation  of 
energy  for  Sc+@C6o  system  for 
LvNMD,  BOMD  and  ADMP 
methods  as  function  of  number 
of  steps. 
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of  endohedral  fullerene  and  transition  metal.  A  prototype  code  has  been  implemented  with 
Gaussian  suite  for  ab  initio  Hartree-Fock  (HF)  and  density  functional  based  tight-binding 
method  (DFTB).  Published  as  paper  10  in  the  archive  below. 

B.  Development  of  Global  Reaction  Route  Mapping  (GRRM)  Method  for  Potential  Energy 
Surface  Exploration.  Efficient  methods  were  developed  for  exploration  of  potential  surfaces  for 
large  molecular  systems  as  well  as  for  finding  the  conical  intersections  and  seams  of  crossing. 
Published  as  Paper  11  and  14. 

C.  R-Matrix  Electron  Scattering  Calculations  for  Diatomic  Molecules 

We  performed  and  published  results  of  R-matrix  calculations  of  differential  cross  sections  for 
low-energy  electron  collisions  with  ground  and  electronically  excited  02  molecule  as  well  as 
integral  and  differential  cross  sections  for  low-energy  electron  impact  excitations  of  N2 
molecule.  Dr.  Tashiro,  who  wrote  the  code  and  performed  calculations,  became  independent 
and  is  continuing  work  along  the  same  line,  and  thus  we  stopped  this  project.  Published  as 
Paper  1  and  2. 

D.  Structure  and  Photochemistry  of  N3CI  and  N3 

Chemistry  of  choloroazide  N3CI  and  hydrogen  azide  N3H  has  attracted  a  substantial 
attention  as  a  potential  source  of  cyclic  N3  species,  the  simplest  nitrogen  high  energy  material, 
as  well  as  in  connection  to  the  NCI-driven  iodine  chemical  laser.  Although  we  studied  the  N3CI. 
Previously  we  theoretically  studied  excited  state  dynamics  of  N3H  as  well  as  the  ground  state 
potential  surface  and  dynamics  of  N3.  In  the  present  grant,  we  focused  on  the  N-CI  bond  energy 
and  photodissociation  dynamics  of  N3CI,  as  well  as  potential  energy  surfaces  of  excited  states 
of  N3. 

i.  Dissociative  photoionization  of  N3CI  and  the  N-CI  bond  energy  in  N3CI.  This  is  a 
collaborative  research  with  Alec  Wodtke’s  high  resolution  synchrotron-radiation-based 
photoionization  mass  spectrometry  experiment.  Using  the  theoretically  calculated  binding 
energy  of  CIN3+  (0.2  eV),  the  adiabatic  ionization  potential  of  CIN3  =  9.97±  0.02  eV  was 
estimated.  The  measured  appearance  potentials  for  N3+,  Cl+,  and  N+  provide  three  independent 
determinations  of  the  Cl-N  bond  energy  in  CIN3,  leading  to  a  new  value  of  the  N-CI  bond  energy 
in  CIN3,  D0(CI-N3)=1 .86  ±  0.05  eV,  0.3  eV  lower  than  previously  reported  values.  The  high  level 
ab  initio  (CCSD(T))  electronic  structure 
calculations  extrapolated  to  the  complete 
basis  set  limit  yield  a  value:  D0(CI- 
N3)=1 .87  eV.  This  work  has  been 
published  Paper  3. 

ii.  Possible  pathway  for  formation  of 
cyclic  N3  in  157nm  photodissociation  of 
N3CI.  The  photodissociation  dynamics  of 
chlorine  azide  (CIN3)  at  the  157  nm  region 
was  studied  theoretically  using  the 
multireference  configuration  interaction 
(MRCI)  method  and  the  complete  active 
space  (CASSCF)  direct  dynamics  method. 

Based  on  the  MRCI  calculations,  the 
excitation  at  the  157nm  region  was 
assigned  to  the  ^A^Sr)  x  ^'(So) 
transition.  A  likely  pathway  for  the 
formation  of  cyclic  N3  after  this  transition 
was  identified  by  direct  CASSCF  dynamics 
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Figure  2.  A  cut  of  potential  surfaces  of  low-lying  N3 
doublet  states  in  C2v  symmetry  as  a  function  of 
ZNNN  at  a  fixed  R(N-N)  of  1 .41  A.  The  large 
black  dot  indicates  the  beginning  of  the  N3 
trajectory,  and  the  large  grey  dots  the  conical 
intersections  encountered.  [sa(4)- 
CASSCF(1 5e/1 2o)/aug-cc-pVTZ  level]. 
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as  follows  (Figure  2).  CIN3  excited  to  ^A^Sy)  dissociates  after  about  40  fs  to  excited  N3(22A’, 
with  about  44  kcal/mol  internal  energy)  +  CI(2P).  This  vibrationally  hot  N3(22A)  goes  diabatically 
through  a  conical  intersection  with  N3(72A’)  at  44  fs  onto  72A’.  At  19  fs  later  and  repeatedly  after 
every  55  fs,  N3(72A’)  crosses  and  trickles  down  via  Coriolis  coupling  to  N3(22A”/2B1)  state,  which 
has  a  potential  minimum  at  the  cyclic-N3  structure.  Some  fraction  of  N3(2^A”/2B1)  produced  will 
survive  dissociation  and  will  be  found  as  the  cyclic-N3,  and  some  other  fraction  will  eventually 
dissociate  to  N(2D)  +  N2  over  a  high  barrier  found  previously.  The  paper  has  been  published  as 
Paper  4. 

iii.  Structures  and  energies  of  low-lying  doublet  excited  states  of  N3  from  accurate 
configuration  interaction  calculations.  The  low-lying  doublet  excited  states  of  the  azide 
radical  (N3)  have  been  studied  at  a  highly-level  multireference  ab  initio  theory  including  basis 
sets  up  to  augmented  quadruple-^  quality.  A  full  hypersurface  scan  under  C2v  restrictions  for 
five  low-lying  bent  N3  states  (2A2,  2A1;  2Bi ,  and  two  2B2)  revealed  a  highly  complex  potential 
surface  manifold  with  many  stationary  points,  conical  intersections  and  multiple  surface 
crossings,  all  of  which  have  been  characterized.  The  behavior  of  these  states  is  discussed, 
especially  as  a  function  of  the  NNN  angle.  At  least  two  new  low-lying  pathways  on  the  excited 
surfaces  leading  from  the  linear  to  the  cyclic-N3  region  were  found,  both  involving  the 
components  of  the  degenerate  excited  2nu  state  of  linear  N3.  The  paper  was  published  as  Paper 
9. 

v.  Analytical  potential  energy  surfaces  for  N3  low-lying  doublet  states.  Adiabatic 
potential  energy  surfaces  (PESs)  for  five  low 
lying  doublet  states  (three  2A’  states  and  two  2A” 
states)  of  N3  are  constructed  based  on  1504 
single  point  calculations  at  the  MRCISD(Q)  level 
using  aug-cc-pVTZ  basis  set.  A  new  strategy  is 
adopted  to  obtain  the  final  PESs  by  combining 
global  fits  of  individual  adiabatic  PESs  and  local 
simultaneous  fits  of  two  adiabatic  PESs  in 
several  conical  intersection  regions  with 
switching  functions.  See  Figure  3,  for  example, 
for  local,  global  and  final  PESs  in  the  vicinity  of 
conical  intersection  1.  These  global  fits  employ 
basis  functions  that  satisfy  permutational 
invariance  with  respect  to  like  nuclei,  and  have 
rms  errors  around  2-3  kcal/mol.  The  special 
local  two-state  fits  are  performed  at  the  cyclic, 
bent  and  linear  N3  conical  intersection  regions  to 
take  account  of  intrinsic  square  root  behavior  of 
the  potentials  and  thus  to  improve  the  quality  of 
fitting.  Stationary  points  as  well  as  minima  on 
the  conical  intersections  and  seams  of  crossing 
are  located  on  these  PESs  and  compared  with 
ab  initio  results:  the  agreement  is  satisfactory  in 
most  cases.  In  addition  to  the  construction  of  adiabatic  PESs,  diabatization  is  performed  for  the 
1  2 A'  and  2  2A’  states  around  their  conical  intersection  at  the  N3  bent  region;  these  two  diabatic 
PESs  and  diabatic  coupling  potential  have  been  constructed  and  reported.  This  work  has  been 
published  as  Paper  8. 

E.  Potential  Energy  Surfaces  and  Dynamics  of  Photodissociation  of  H2CO 


(0  (I) 

Figure  3.  Colnl  potential  energy  surfaces. 
(a),(d):  surface  from  global  fit;  (b),(e):  surface 
from  local  fit;  (c),(f):  final  surface,  rl  =  r2  for 
(a),  (b),  (c);  <NNN  =  60s  for  (d),  (e),  (f). 
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Formaldehyde  has  provided  many  opportunities  to  both  experimentalists  and  theoreticians 
to  study  the  mechanism  of  photodissociation  as  a  simplest  and  prototypical  molecule.  However, 
the  mechanism,  in  particular,  how  the  molecule  excited  by  radiation  to  the  first  singlet  excited 
state  (Si)  eventually  comes  down  to  the  ground  state  (S0)  to  dissociate  mainly  to  H2  +  CO.  A 
burst  of  theoretical  activities  has  taken  place  in  the  last  few  years,  and  we  have  been  one  of  the 
contributors  to  this  exciting  story. 

i.  Photodissociation  dynamics  of  formaldehyde  initiated  at  the  IVSo  minimum  energy 
crossing  configuration.  The  photodissociation  dynamics  of  H2CO  is  known  to  involve  Si,  Ti 
and  S0.  Recent  quasiclassical  trajectory  (QCT)  calculations  in  conjunction  with  experiment  have 
identified  a  “roaming”  H-atom  pathway  to  the  molecular  products,  H2+CO.  These  calculations 
were  initiated  at  the  global  minimum  (GM)  of  S0,  which  is  where  the  initial  wavefunction  is 
located.  The  “roaming”  mechanism  is  not  seen  if  trajectories  are  initiated  from  the  molecular 
transition  state  saddle  point  (SP).  In  this  work  we  identify  the  minimum  energy-crossing 
configuration  and  energy  of  the  Ti/S0  potentials  as  a  step  towards  studying  the  multi-surface 
nature  of  the  photodissociation.  QCT  calculations  are  initiated  at  this  configuration  on  a  revised 
potential  energy  surface  and  the  results  are  compared  to  those  initiated,  as  previously,  from  the 

50  GM  as  well  as  the  S0  SP.  This  work  has  been  published  as  Paper  5. 

ii.  Automated  global  mapping  of  minimal  energy  points  on  seams  of  crossing  by  the 
global  reaction  route  mapping  (GRRM)  method:  A  case  study  on  H2CO.  The  automated 
global  reaction  route  mapping  (GRRM)  method  using  the  anharmonic  downward  distortion 
following  (ADDF)  procedure  has  previously  been  developed  and  applied  to  single  potential 
energy  surfaces  (PESs).  Dr.  Satoshi  Maeda  who  developed  the  method  recently  joined  our 
group  as  Visiting  Research  Fellow.  We  have  extended  method  to  perform  automated  global 
mapping  for  minimal  energy  points  on  seams 
of  crossing  (MSX  structures).  Here  the  ADDF 
is  applied  to  a  penalty  function  based  on  two 
PESs  of  different  electronic  states.  The 
present  approach  is  effective  not  only  for 
crossing  seams  between  states  with  different 
symmetry  but  also  for  conical  intersections  for 
states  with  the  same  symmetry.  As  shown  in 
Figure  4,  many  new  MSX  structures  were 
discovered  on  the  So/T)  and  S^/T^  crossing 
seams  and  the  S0/Si  conical  intersections  of 
H2CO  by  automated  global  mapping  using  the 
ADD-following  method.  A  possible  pathway  for 
dissociation  of  formaldehyde  excited  to  Si  at 
low  energy  is  discussed.  This  work  was 
published  as  Paper  7. 

iii.  Photochemical  reactions  of  the  low- 
lying  excited  states  of  formaldehyde.  T,/S0 
intersystem  crossings,  characteristics  of  the 

51  and  T,  potential  energy  surfaces,  and  a 
global  Ti  potential  energy  surface.  Accurate 
ab  initio  calculations  using  the  multireference 
configuration  interaction  method  have  been 
performed  to  characterize  the  potential  energy 
surfaces  (PESs)  of  low-lying  excited  states  (Si 
and  Ti)  of  formaldehyde  (H2CO)  and 
hydroxymethylene  (HCOH)  with  emphasis  on  their  isomerization,  dissociation  and  the  possible 
role  of  the  7)  state  in  the  non-adiabatic  photodissociation  of  H2CO,  as  shown  for  T,  in  Figure  6. 


Si/Ti-MSXI 

-»430.9  kJ  mol-'  V7VMSX1 
0.0  kJ  mol'1 


H2  +  GO 

218.9  kJ  mol  1 

planar  eA-HCOH  planar  /ram-HCOH 

Figure  4.  Illustration  of  the  potential  energy 
profiles  for  inside  the  potential  wells  of 
formaldehyde  (A)  and  hydroxymethylene 
(B)  based  on  the  CASPT2  theory. 
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Two  regions  on  the  7)  PES  are  found  to  contribute  to  the  non-adiabatic  transition  to  the  ground 
(S0)  state.  Three  minima  on  the  seam  of  crossing  (MSXs),  80-85  kcal/mol  (above  the  S0  global 
minimum),  are  located  in  the  HCOH  region;  they  however  are  blocked  by  a  high-energy 
isomerization  transition  state  at  -107  kcal/mol.  The  other  MSX  discovered  in  the  H2CO  region  is 
reachable  with  energy  <  91  kcal/mol  and  strong  spin-orbit  interaction;  this  may  be  a  more 
important  pathway  for  the  7)  to  S0  transition.  A  full-dimensional  PES  is  generated  for  the  7) 
state,  fitted  by  a  weighted  least  squares  method  employing  a  many-body  expansion  in  which 
each  term  is  a  function  of  the  internuclear  distances  and  is  invariant  under  permutations  of  like 
atoms.  The  single  global  function  covers  the  formaldehyde  and  the  hydroxymethylene  regions 
as  well  as  dissociation  pathways.  The  high  quality  of  the  fitted  PES  is  demonstrated  by  the 
small  root-mean-square  fitting  error  of  119  cm-1  and  the  close  agreement  between  the  critical 
points  from  ab  initio  calculations  and  from  the  fitted  PES.  This  work  was  published  as  Paper  8. 

F.  The  Mechanism  of  ion-Molecule  Reaction:  POCI3'  +  02 

In  the  present  experimental-theoretical  joint  efforts,  the  oxidation  of  the 
trichlorooxyphosphorus  anion  (POCI3“) : 

POCI3“(  X2A' )  +  02(  X3Xg )  ->  P02CI2-(  XM' )  +  CIO(  X2n )  +  45.4  kcal/mol 

which  takes  place  in  combustion  flames  has  been  examined  experimentally  in  a  variety  of 
temperatures  and  theoretically  via  ab  initio  and  density  functional  methods.  The  reaction  was 
studied  experimentally  by  the  team  of  Dr.  Al  Viggiano  using  a  turbulent  ion  flow  tube  (TIFT)  and 
kinetics  was  measured  between  300  and  626  K  estimating  an  overall  reaction  barrier  and 
reaction  exothermicity  of  1.23  and  45.4  kcal/mol,  respectively.  While  a  plethora  of  experimental 
data  has  been  accumulated  in  the  last  few  years  about  the  reactivity  of  various  oxyphosphorus 
compounds,  there  are  virtually  no  detailed 
studies  available  that  look  into  their  reactivity 
from  a  theoretical  standpoint.  We  have 
performed  DFT  and  ab  initio  calculations  of 
the  potential  energy  surface  involved  in  the 
reaction.  As  shown  in  Figure  2(A),  we  found 
at  the  UB3LYP/6-31+G*  level  that  the 
reaction  proceeds  in  the  doublet  state  with  a 
multi-step  mechanism,  in  which  reactants  at 
first  form  a  pre-reaction  ion-molecule 
complex  intermediate  (Rcompl)  and  go  over 
a  barrier  (TSO)  to  reach  a  five-coordinate 
trigonal  bipyramidal  (TBP)  intermediate 
[CI3P(0)(00)]‘  (LMO),  which  through  a  very 
skewed  substitution  (breaking  a  Cl-P  bond 
and  forming  a  CI-0  bond)  transition  state 
TS1  forms  a  chloroperoxy  intermediate 
[CI2P(0)(00CI)]"  (LM1).  Finally  the  00  bond  is  broken  in  LM1  via  TS2  to  form  the  product 
complex  [CI2P(0)0]"(CI0)  (Pcompl),  which  finally  dissociates  to  the  products  [CI2P(0)0]"  +  CIO. 
Among  various  possible  positions  of  ligands  in  the  five-coordinate  TBP  complex,  only  two  forms 
and  thus  two  reaction  pathways  have  been  found  to  exist.  The  equatorial  pathway  involves 
attack  of  02  cis  to  the  PO  bond  in  [CI3PO]"  giving  the  TBP  intermediate  [CI3P(0)(00)j"  (LMO)  in 
which  the  O  ligand  is  in  the  equatorial  plane  (see  Fig.  2(A)),  and  the  axial  pathway  has  trans  02 
attack  and  the  TBP  intermediate  with  the  O  ligand  (as  well  as  the  00  ligand)  in  the  axial 
position  (not  shown).  The  rate-determining  barrier  is  at  TS1  and  is  about  the  same,  12.2  and 
12.0  kcal/mol  for  equatorial  and  axial  pathways,  respectively,  at  the  UB3LYP/6-31+G*  level,  and 
two  pathways  should  compete  against  each  other.  DFT  is  knows  to  be  unreliable  in  describe 
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charge  delocalization.  Higher  level  single  point  calculations  (shown  in  Figure  2(B))  lower  the 
barriers  substantially  and  give  the  barrier  height  for  the  equatorial  pathway  to  be  3.2  and  2.5 
kcal/mol  at  TSO  and  TS1,  respectively,  at  the  UCCSD(T)//aV(D+d)Z  level,  in  excellent 
agreement  with  the  experiment.  Radical-radical  reactions  usually  have  no  or  very  small  barrier. 
The  (small)  barrier  in  the  present  reaction  seems  to  be  caused  by  the  energy  required  to  deform 
the  nearly  tetrahedral  reactant  structure  to  a  trigonal  bipyramid  structure  at  TSO  as  well  as  small 
electron  affinity  of  02  (0.45eV).  The  work  was  published  as  Paper  12. 

G.  Reactions  of  VX  Agent 

(i)  Proton  affinity  and  fluoride  affinity  of  the  VX  toxic  nerve  agent.  Protonation  and 
fluorination  were  found  to  give  various  unique  structures  including  dissociated  species  and  ion 
molecule  complexes. 

(ii)  Chemical  ionization  detection  for  monitoring  of  VX  simulants  with  CIMS.  An 

experimental  and  theoretical  study  was  published  as  Paper  8. 
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